Far-off-resonant wave interaction in one-dimensional photonic crystals with quadratic 

nonlinearity 

Johannes-Geert Hagmann *, Lasha Tkeshelashvili, and Kurt Busch 
Institut fur Theoretische Festkorperphysik, Universitat Karlsruhe (TH), 76128 Karlsruhe, Germany 



Guido Schneider 
Institut fur Analysis, Dynamik und Modellierung, 
Universitat Stuttgart, Pfaffenwaldring 57, 70569 Stuttgart, Germany 
(Dated: February 6, 2008) 

We extend a recently developed Hamiltonian formalism for nonlinear wave interaction processes 
in spatially periodic dielectric structures to the far-off-resonant regime, and investigate numerically 
j the three-wave resonance conditions in a one-dimensional optical medium with nonlinearity. In 

particular, we demonstrate that the cascading of nonresonant wave interaction processes generates an 
effective x' 3 ' nonlinear response in these systems. We obtain the corresponding coupling coefficients 
' through appropriate normal form transformations that formally lead to the Zakharov equation for 

pP , spatially periodic optical media. 
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I. INTRODUCTION 



Oh| During the last two decades, the emergence of periodically structured optical materials - commonly called photonic 
crystals (PCs) - has lead to substantial progress in the science and technology of optics and photonics. For a recent 
review, we refer to Ref. [I]. It has soon been realized that PCs with nonlinear constituent materials facilitate 
the realization of novel effects and functionalities that are very hard to realize with other systems. However, a 
mathematically consistent description of wave-interaction processes in PCs is complicated due to a number of formal 



and as yet unsolved problems relating to the complex structure of the equations of motion. Recently, based on a 
Hamiltonian formulation of Maxwell's equations, an effective field theory for nonlinear PCs has been developed by 
Volkov and Sipe [2| • This approach allows a systematic study of nonlinear effects beyond the commonly used multiple- 
scales analysis Q . In particular, in Ref. [|J , the formalism has been applied to the problem of resonant interactions of 
wave packets in PCs with quadratic nonlinear response. The purpose of this article is to extend this formalism to the 
far-off-resonant, i.e. nonresonant, case in one-dimensional systems. In particular, we reduce the nonlinear equations 
' of motions of corresponding nonlinear PCs to the simplest form possible. This reduction is facilitated through the 
, systematic elimination of nonresonant interaction terms from the classical Hamiltonian via appropriate normal form 
transformations. This procedure results in the emergence of an effective third-order response of the system. We would 
, like to note that similar techniques have been successfully employed in the analysis for magnetic systems [J] as well 



as for ideal fluid hydrodynamic waves [f|. In the latter case, the result is the so-called Zakharov equation for surface 
OO ! water waves, which has recently been studied experimentally @. 

However, in contrast to the case of spin or water waves, where the underlying linear dispersion relations prohibit 
nontrivial three- wave resonances, such resonances do appear in PCs for certain sets of interacting wave numbers. This 
leads to nonremovable singular terms in the normalized Hamiltonian. Consequently, additional constraints on the 
wave packets are required in order to mathematically justify the validity of the simplified equations for finite times. In 
realistic physical systems, where losses due to fabricational tolerances and/or residual material absorption, or simply 
the finite sample size provide natural upper bounds on the time scales, these constraints may be less severe. 

The article is organized as follows. In Sec. QT] we extend the Hamiltonian formalism of Volkov and Sipe [|J for 
one-dimensional photonic crystals with \^ nonlinearity to the nonresonant case. In Sec. IIII[ we analyze numerically 
the lowest-order resonance conditions and compare the results to the water wave problem for which the spectrum 
is known analytically. In Sec. IIV1 the normal form reduction of the Hamiltonian is performed and the Zakharov 
equation for nonlinear PCs is formally derived for the case when three-wave interaction processes are suppressed. In 
particular, the expressions for the effective third-order nonlinear response are derived and compared with the case of 
homogeneous dispersive x*- 2 ) materials. We conclude with a discussion of the limitations of our approach, and provide 
an outlook on possible future work in Sec. |Vj 
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II. CLASSICAL HAMILTONIAN FORMULATION 



In this section we derive a classical Hamiltonian formulation for nonlinear PCs. The derivation closely follows Ref. 
0, however, we emphasize that, as compared to Ref. @, we have chosen a different set of fundamental fields for the 
decomposition into a complete set of eigenf unctions. 

We consider infinite nonlinear PCs with nonmagnetic, dispersion- and lossless constituent materials. Since the linear 
dispersion provided by photonic crystals' microstructure is typically much larger than the material dispersion of its 
constituent materials, this is a very good approximation. If material dispersion dominates - such as would be the case 
for quasiphase matched gratings where only the nonlinear properties vary periodically in space - one would have to 
resort to a slowly varying envelope approximation Q. Assuming further that the material response is local in space 
and time, we can express the polarization P(x, t) that connects the dielectric displacement D(x, t) = eoE(x, i)+P(x, t) 
with the electric field E(x, t) as 

A(x,i) = e ^(x,<) + x!) ) (x)^(x,<)+ X ^(x)£; 3 (x,i)£; fc (x,i) + ... . (I) 

Here, eo denotes the free-space dielectric permittivity. The susceptibility tensors x 01 rank n + I are symmetric 
under permutations of their indices and - here as well as in the remainder of the article - summation over repeated 
indices is implied. Instead of the first-order susceptibility x^\ we employ the dielectric tensor e(x). In addition, in 
the following, we will restrict ourselves to the case of isotropic linear response, so that we have %W(x) = eo ( e ( x ) — 1). 

Throughout this paper, we identify the canonical momentum with the displacement field D rather than with the 
electric field E. Therefore, it is desirable to expand the electric field in terms of the displacement field according to 

£?<(x,t) -e - 1 e(x)- 1 A(x,i)+rgl(x)Z?,(x,<)C fe (x,i) + ... . (2) 

For the second-order susceptibility, the relations between x^ and T^ 2 ) reads 

y (2) fxl 

^ j _ (eoe(x)) 3 ' [ } 

and the relations for the higher-order susceptibilities can be found by recursion 0]. In PCs, the material properties 
are periodic functions with respect to the set of lattice vectors R, i.e., we have e(R + x) = e(x), xtu l \ (R + x) = 

X&1(*), ifc^CR + x) = lfc.^(x), etc. 

Within the Dzyaloshinski or generalized Coulomb gauge where the scalar potential <p(x, t) = 0, we may express 
both, the electric field E and the magnetic field B, through the vector potential A according to 

E(x,f) = -d t A{x,t), (4) 
B(x,t) = V x A(x,<). (5) 

For this system, the Hamiltonian that yields Maxwell's equations in the absence of external sources and currents reads 

H(A,n) = ~ y"d 3 x^-^-yA(x,t)A(x,t)+^o 1 -B l (x,i) J B l (x,<) 

+ 1 J d 3 xTff k (x.)D i (x,t)D j (x,t)D k (x,t) + ... 
= H 2 +H z + ... . (6) 

In the former expression, II = D denotes the canonical momentum associated with the vector potential A 0]. 
Consequently, the equal-time Poisson brackets are 



{ ^(x,t),n j( x',t)} An = /dV'f##^g#4 



8Ai{x,t) <m.,-(x',i) 



t) 5Il k (x",t)6A k (x",t), 
= %<5(x-x'), (7) 

where S/5Ai(x, t) and S/SHi(x,t) denote functional derivatives (see Ref. 0] for a definition and properties). The - 
on first sight somewhat superfluous - indices A and II on the left-hand side allow us to keep track of the fields, once 
a transformation is applied to the Hamiltonian ©. To decompose this Hamiltonian into symmetry-adapted basis 
functions (modes), it is useful to consider the solutions of the linearized equations of motion for periodic systems. 
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A. Bloch function decomposition 



In the linear limit, D(x, t) = eo e ( x )E(x, t), and the Maxwell equations reduce to two decoupled partial differential 
equations for the fields A(x, t) and D(x, t). Since we are interested in stationary solutions (modes), we employ the 
time-harmonic ansatz A(x, t) = A nk (x) exp(— iw^t), D(x, t) = D n k(x) exp(— iWnkt), and obtain from the Maxwell 
equations 
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■V x (V x A nk (x)) - -f A nk (x 



D Ilk (x)\\ _ LO 



,2 



v H Vx b£nJ = ^ Dnk(x) - (8) 

In the above ansatz, we already made use of the Floquct-Bloch theorem that - in the case of lattice periodic functions 
e(x) - applies to Eqs. © and, therefore, introduced the Bloch functions A nk (x) and D„i c (x). These Bloch functions 
are labeled through the discrete band index n > and the continuous wave vector k G Jl^ that lies in the first Brillouin 
zone fif,. This labeling is commonly referred to as the reduced zone scheme [llj . In general, the Bloch functions and 
their associated eigenvalues u nk have to be determined numerically through one of the well-documented techniques 
1]. For future reference, we would like to note that the Bloch functions have the forms of modulated plane waves 

D„ k (x) = d„ k (x)exp(ikx), (9) 
A nk (x) = a„ k (x)exp(zkx), (10) 

where a nk (x) and d„k(x) are lattice-periodic functions. For convenience, the Bloch functions D„k, A„k and the 
eigenvalues are chosen to be periodic in reciprocal lattice vectors G, i.e. D„(k+G) = D nk for all reciprocal lattice 
vectors G (T^ |. It is also useful to note that the time-reversal symmetry of our problem implies cj„k = Wn-k which, 
in turn, implies that together with the solutions (D„ kl V x A„ k ), the corresponding pair (D*_ k , — V x A*_ k ) also 
fulfills the Maxwell equations for the same frequency. Therefore, we choose D nk = D*_ k which guarantees that the 
displacement field D(x, t) is real. The Bloch functions constitute a complete orthonormal set of basis functions. Since 
Eq. implies 

D„ k (x) = ie e(x)w„ k A nk (x), (11) 

we find the the orthogonality relations [ici| 

| Af e(x)w„ k A„ k (x)A*, k ,(x) = i<W£(k-k'), (12) 
d 3 x(Vx A nk (x))-(Vx A*, k ,(x)) = W^±6 nn ,6(k-k>), (13) 



for k, k' € tts if n ^ 0. For n = 0, k = 0, we have w rlk = and the normalization has to be performed separately. 
This issue will be addressed below. 

For the subsequent formal developments, we find it more transparent to represent the Bloch functions in the so- 
called extended zone scheme [Tl| . where the wave vector k varies over the entire reciprocal space. Formally, this is 
facilitated by defining the band index to be a function of the wave vector, i.e., to introduce 

A k (x) = A„ (k)k (x), (14) 

and an analogous definition for D k and cD k - Whether such a mapping does or does not exist depends on the spectrum for 
the selected system, in particular on the occurrence of degeneracies. For one-dimensional PCs and wave propagation 
along the axis of periodicity, wave vector and Bloch functions reduce to scalars. As a result, the unfolding of 
the band structure is realized through n{k) := [\k\d/n\ e N, where |_-J denotes the Gauss bracket defined via 
< \k\ — {■K/d)n(k) < (ir/d) (d denotes the lattice constant of the PC). In two- or three-dimensional systems, this 
transformation between reduced and extended zone scheme is not straightforward but can, in principle, be carried out 
numerically. Therefore, we in the following restrict ourselves to the case of strictly one-dimensional PCs and assume 
that the system is periodic in the x direction. To simplify the notation, we assume in the following an isotropic 
nonlinear response of the constituent materials so that the relevant fields (A, D, and E) are all polarized along a 
certain transverse direction which we define as the y direction. As a result, we may suppress the vectorial character 
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of all fields. In such a one-dimensional PCs, Eqs. (| 1 1 1) - (fT3f read as 

D k (x) = ie e(x)u> k A k (x), (15) 
dx eoe(x) iD k A k (x)A* k ,(x) = -S(k-k'), (16) 

dx d x A k {x) d x A* k ,(x) = ^5{k- k'), (17) 
for k,k' e R\{0}. 

We now return to the above-mentioned normalization problem related to u) k= a — for which the eigenvalue equation 
reads as 

d xx D (x) = 0. (18) 

Owing to the Bloch-Floquet condition and our choice of D k = D*_ k , we obtain Do(x) = a € R. On the other hand, 
D k must be a continuous function of k so that the scaling D k oc y/uk requires us to define 

D k = Q = 0, (19) 
lira Dn oc lim JJi. (20) 

Therefore, an absolute normalization in the sense of the orthogonality relations is impossible, yet our definition is 
justified from the physical point of view [HI]. In fact, physically the mode Q k =0 = corresponds to a simple translation 
of our infinite PC and this clearly requires an infinite amount of energy. Consequently, this mode resists an absolute 
normalization. The above definition ensures that this is the case. Clearly, this discussion may be directly transferred 
to the case of two- and three-dimensional PCs. Finally, we expand the physical fields D(x,t) and A(x,t) of the 
nonlinear system into the orthogonal basis of the linearized problem 



D(x,t) = / dkD k (x)a k (t) + c.c, (21) 
Jm. 

A(x,t) = / dkA k (x)a k (t) + c.c, (22) 



where the time-dependent amplitudes a k (t) and a^(t) are the weights associated with each Bloch mode. In terms of 
the physical fields, the new set of variables (k ^ 0) are expressed as 

a k (t) = [dx f_L D *k(x)D(*,t) + eQe{x) Qk A* k (x)A(x, t)) . (23) 

In the derivation of Eq. (|23p . the identity (D k , d x A k ) — (D*_ k ,—d x A k _ k ) has been employed and it can be verified 
that Eq. (|23[) represents a canonical transformation. Consequently, for arbitrary functions F, G in phase space, 



{F,G} DtA = {F > G}- a - a *. (24) 
holds. In terms of the new variables, the equal-time Poisson brackets now have the explicit form 

^ , ^ aa J R '(<5afe(5aJ Sa k Sal) ' ^ ^ 

so that 

{a k ,a* k ,}a,a* = iS(k-k'), (26) 
{a k ,a k >}a,a* = 0, (27) 

is immediate. 



B. Hamiltonian in the new variables 



The results of the previous section may now be utilized to rewrite the original Hamiltonian ^ in terms of the new 
field variables a k and a£. The linear part, H2, is a quadratic form 



H 2 = J dk oj k a k a* k , 



(28) 
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while the lowest-order nonlinear term, H 3 , couples three Bloch modes and reads as 

H 3 = 



/ dk 123 V^l 2k;j a kl a k2 a k3 S(G ~k x -k 2 - fc 3 ) 

J G 

f / dk 123 V k ^l 2k3 a* ki a k2 a k3 ^2d(G + k 1 -k 2 - k 3 ) 



+ c.c. (29) 

= H p + H m + c.c. 

In the former expression, the integration is carried out over all three momenta, i.e., dki 23 = dk x dk 2 dk 3 and the sum 
over the entire set of reciprocal lattice vectors G £ {G el: G = m(2w/d) , m £ Z} reflects the fact that, in 
the extended zone description, in principle all bands are coupled. In the reduced zone scheme, this would include 
both the so-called Umklapp processes, i.e., when the sum of the coupled wave vectors leaves the Brillouin zone as 
well as the coupling between different bands without Umklapp processes. Indeed, this complication of classifying and 
keeping track of the various scattering processes in the reduced zone scheme has lead us to adopt the (perhaps less 
familiar) extended zone scheme. The corresponding coupling coefficients V k ^ k2k3 and V^ 2 ^ are derived from the 
Bloch functions as given in Eqs. © and (TIT)]) (transferred to the extended zone notation) together with the identity 

^exp(-ifci?) = ^^5(G + fc), (30) 

R G 

that allows one to reduce the integration over all space into an integration over the PC's unit cell V cs n 



KLk s = TT I dxT^(x)D kl (x)D k2 (x)D k3 (x) , (31) 

-Veen 



« _ £ZL / ^ r(2) 



vSlks = %[ Vc j v dxT^(x)Dt 1 (x)D k2 (x)D k3 (x). (32) 

In this expression, we have introduced the cyclic permutation operator V c that acts on any function F\ 23 with indices 

123 according to V C F\ 23 — F\ 23 + F 2 \ 3 + F 32 i. Physically, the terms associated with V^^ afc correspond to the 

( 2 \ 

simultaneous creation or annihilation of three waves, while the terms associated with VL L u describe wave mixing 

' KiK 2 K 3 o 

processes where a single wave is converted into two other waves. Clearly, a process associated to the matrix element 
^feifefe does not conserve energy. In order to be efficient, a coupling mechanism to another physical system is required. 
Therefore, for conservative optical systems these processes are strongly suppressed. 

The explicit form of the higher-order nonlinear contributions to the Hamiltonian, Hj, j > 3, follow in strict analogy 
to the derivation of H 3 . For instance, the next nonlinear term Hi contains both energy-suppressed and energetically 
allowed processes. Among all of those, the relevant part of the Hamiltonian corresponding to the third-order nonlinear 
susceptibility is 

Hl 3) = f dk 12U W^ k3kA a kl a k2 at 3 a* ki ]T S(G + k x + k 2 -k 3 ~ k 4 ), (33) 

J G 

where processes in which two waves are scattered are described via the coupling coefficient 

W SU 3ki = I dx r {3 Hx)DkAx,t)D k2 (x,t)Dl(x,t)Dl(x,t). (34) 

za Jv CBll 

Here, the permutation operator Vq acts on a function with indices -F1234 according to VqFi 23 4 = -F1234 + -F1324 + 

^3214 + -Fl432 + -F4231 + -^4321- 

The equations of motion for the canonical variables are the Hamilton equations 

a k (t) = {d k (t),H}- ara <. (35) 

However, it is often preferable to work in a frame that "rotates" with the frequency of the free carrier wave, i.e., to 
employ the slow variables a k (t) — d k (t) exp(—iCu k t) . In this case the equations of motions read 

d t a k = l /*23^ fc 2 2 ) fe3 afc 2 afc3e l( " fc2+ "' fc3 ^ )t ^^G + fc-fc2-fc3) (36) 

J G 

+{a k {t),H-Hf ) }- ara . (37) 
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The terms on the right-hand side of the equation of motion give an oscillatory contribution to the mode unless 
certain resonance conditions are met. For the one-dimensional photonic crystals, these conditions are analyzed and 
compared to a related problem in ideal hydrodynamics in the following section. 



III. THREE- WAVE INTERACTION 



In many problems of interest in optics, the nonlinearities that are involved in the system dynamics are quite small. 
This suggests that in the equations of motion the nonlinear interaction terms can be considered as a perturbation to 
the linear part. The resonance conditions for nonlinear wave mixing processes are determined by the linear dispersion 
relation u>k ■ A profound example of such a resonant mixing process in optics is the second harmonic generation, where 
the frequency of an incident beam is doubled through a quadratic, i.e., x*- 2 -* -nonlinear response of the system. The 
requirement for efficient conversion is the conservation of energy and the momentum in the process, lu + ui = 2uj and 
k + k = 2k, respectively. These equations are examples of the more general class of three- wave resonance conditions 

siuj(ki) + s 2 Lu(k 2 ) + s 3 Lu(k 3 ) = (38) 
siki + s 2 k 2 + s 3 k 3 , = (39) 

where Si = ±1. These conditions are the requirements for an efficient energy transfer between three waves (modes) 
of wave numbers ki,k 2 ,k 3 . For PCs, the second condition must be replaced by the conservation law for crystal 
momentum, i.e., reciprocal lattice vectors can be added to any of the momenta. For the moment, we leave this point 
aside but will return to it at a later stage of the discussion. 

A detailed discussion of selection rules for nonlinear PCs can be found in Ref. [14]]. We would like to note that in the 
case of non- negligible third-order response, i.e., x nonlinearities, similar conditions must hold for four- wave mixing 
processes. At first sight, it appears difficult to solve Eqs. (|3"9"j) even if the linear spectrum lu^ is known analytically 
such as in the case of surface waves in an ideal fluid. On the other hand, for fixed Si, Eqs. (|39|) contains six variables 
but only three of them can be chosen independently, and the two equations remove one degree of freedom. In this 
reduced parameters space and for any given tOk , one can easily carry out the analysis numerically and visualize the 
resonance conditions. 

Upon combining the expressions in Eqs. (|39[) . we obtain 

Lu(-sis 2 k 2 - sis 3 k 3 ) + sis 2 uj(k 2 ) + sis 3 uj(k 3 ) = (40) 

Equation ([40]) is invariant under the simultaneous inversion of the signs of all Si , so that only four out of eight cases 
have to be considered. Time reversal symmetry, i.e., oj(k) = uj(—k) > 0, and the fact that in one-dimensional PCs 
the spectrum is nondegenerate, lead to further reduction. In fact, a semi-positive spectrum excludes the case of three 
identical signs and the time-reversal symmetry maps two other cases onto each other up to an exchange of indices. 
Consequently, in order to find nontrivial solutions to Eq. P0|) . we are left with only two distinct cases. These can 
readily be visualized by defining the two functions 

lu u \ ^( fc 2 - k 3 )-uj(k 2 ) 

ai(k 2 ,k 3 ) = — h 1 (41) 

iu{k 3 ) 

lu u \ u(k 2 + k 3 ) -uj{k 2 ) , 
a 2 (k 2 ,k 3 ) = — h 1 (42) 

-U}(k 3 ) 

for k 2 £ K, k 3 € K \ {fco = 0}, (k 3 = fc only yields a trivial resonance). With these definitions, the resonance 
conditions are fulfilled if and only if a,i{k 2 , k 3 ) = 0. 

The first spectrum tok for which we evaluate Eq. (|42[) is the dispersion relation for gravity-capillary waves on the 
surface of an ideal fluid with infinite depth 



^ fe = ^g\k\+a\k\\ (43) 



where g is the gravitational acceleration and a the capillarity of the fluid. Introducing k' — yog *fc (see also the 
right panel of Fig. [I}, transfers the previous expression into 

Wk^V* = VW\+W (44) 

for evaluation in ai(k 2 , k 3 ), a 2 (k 2 , k 3 ). Notice that both uif.ia*g~i and k' are dimensionless. It is, therefore, possible 
to discuss the position of the resonances as a function of the gravitational and the capillary constant. The left panels 



FIG. 1: Left panel : \ai(k'2, k'^)\ for gravity-capillary waves; Right panel: Corresponding dispersion relation u)y . All variables 
are dimensionless, for notation and normalization see the text. 




FIG. 2: Left panel : \a2(k' 2 ,k' 3 )\ for gravity-capillary waves; Right panel : Close up showing the exact resonance near k' 2 = 0. 
All variables are dimensionless, for notation and normalization see the text. 



of Figs. [T] and [3 reveal the position of resonant or nearly resonant sets of wave vectors k! 2 £ K, k 3 £ K \ {0}, i.e., sets 
for which \ai(k' 2 ,k' 3 )\ < 0.01. Black and grey regions correspond to \ai(k' 2 ,k' 3 )\ <C 0.01 or exactly zero, whereas white 
regions delineate wave vectors for which \ai(k' 2 , k' 3 )\ > 0.01. The symmetry with respect to the origin is a consequence 
of time reversal symmetry which implies that ai(k' 2 ,k' 3 ) — ai(—k' 2 , — k 3 ). Besides the trivial resonances k' 2 = k 3 for 
ai(k' 2 ,k' 3 ) (Fig. [T]) and k 2 = for ai(k' 2 ,k' 3 ) (Fig. [5]), there are other sets of resonant and nearly resonant pairs of 
wave vectors that are located on ribbons. It is instructive to notice the influence of the capillary constant on the 
position of these ribbons by a scaling argument: Selecting fixed (k' 2 , k 3 ) in the set of nearly resonant wave numbers, 
the location of the point moves in the direction of increasing (k 2l k 3 ) with decreasing a such that in the limit a — ► 
these resonances lie in the infinite for (k 2 , k%). For pure gravity waves {p = 0), there are only trivial resonances when 
one of the wave vectors involves the dc component (u>o = 0). The first nontrivial mixing processes occur at higher 
order (four- wave processes). 

Returning to three-wave mixing in PCs, the situation is very different from the previous case. In contrast to water 
waves, there is no explicit expression for lo^ and the discussion is necessarily qualitative. In order to illustrate the 
situation, we assume that the PC consists of alternating layers with equal thickness but different refractive indices, ua 
and riB, respectively. Then, the only parameter that alters the spectrum is the refractive index contrast An = tia — ns- 
To calculate the photonic band structure, Uk, we use a plane- wave expansion method. Of course, any other method for 
obtaining photonic band structures could be used [l[. Except for the case of a homogeneous dielectric, the spectrum 
takes the form of a band structure, as shown in right panel of Fig. [3J An example with an index contrast of An = 2 
is shown in Figs. [3] and [4] and gives important qualitative information about the occurrence of resonances (dark 
regions) and quasi resonances (grey scale). In these figures, the wave vectors are multiplied with d/2ir, consequently, 
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FIG. 3: Left panel : \a 1 (k2, k^)\ for a PC with refractive index contrast An = 2; Right panel: Corresponding photonic band 
structure uik in the extended zone scheme. 
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FIG. 4: Left panel : 02(^2 , ^3) | for a PC with refractive index contrast An = 2; Right panel: Close up of the upper right 
structure of the left panel. 



the borders of successive Brillouin zones appear at half-integer values of the dimensionless variable. As in the case 
of water waves, the existence of uio — leads to trivial resonances in the diagonal of Fig. [3] and a horizontal line in 
Fig. [4] In addition, there are isolated resonant and nearly resonant subsets of wave vectors; their number and relative 
position depends on the refractive index contrast. However, it is hard to quantify this dependence. In the extreme 
case of a homogeneous material, An = 0, the resonance conditions can trivially be fulfilled. For instance, for 02(^2, k^) 
the two squares ki > 0, k% > and &2 < 0, k$ < will become identically zero. With increasing contrast, photonic 
band gaps form at the borders of the first Brillouin zone and the deviation from a linear spectrum is strongest in this 
region. This behavior is reflected by the pattern in the left panel of Fig. 0J As alluded to above, it is important to 
recall that in PCs only crystal momentum is conserved, while the Oi(/c2, fe) in Figs. [3] and 2] have been calculated on 
the basis of momentum conservation. Therefore, for a complete interpretation all points in the plane (&2, fe) have to 
be shifted through the reciprocal lattice vector G = m2ir/d, m G Z, both vertically and horizontally. Consequently, 
there are more resonant zones than those shown in Figs. [3] and [4] that have to be considered. 

Although at this point, we have not yet formally constructed the normal-form Hamiltonian for the PC, we can 
already anticipate an important difficulty in its derivation. The resonance conditions, Eqs. (|40[) . will appear in 
the denominators of the effective interaction coefficients. Whereas singularities resulting from trivial resonances are 
removable, the singular terms associated with nontrivial resonances are pertinent. Therefore, additional constraints 
on the distribution of wave vectors are required in order to ensure the applicability of the normal form transformation. 
Since a wave packet in a nonlinear medium generally acquires new wave vector components as time goes on, these 
constraints will translate into constraints on the time scales for which our theory will be applicable. 
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IV. REDUCTION TO NORMAL FORM 

Normal form transformations represent a particularly useful tool for the study of nonresonant wave interaction 
processes in Hamiltonian systems. In essence, they facilitate the reduction of nonlinear dynamical systems to the sim- 
plest possible form by eliminating inessential terms or degrees of freedom from the original Hamiltonian. For physical 
systems with continuous degrees of freedom, a good illustration of this principle is the study of the aforementioned 
water wave problem which dates back to the early work of Zakharov [H . More recent highlights include the proof of 
the nonintegrability of the water wave problem with the help of normal form transformations by Craig and Worfolk 

HQ. 

In this section, we employ normal form transformations to the case of nonresonant wave interaction in \^ PCs. 
The goal is derive an effective x' 3 ^ response of the system that originates from a cascading of \^ processes. This 
analysis shows how the engineering of the PC structure could, in principle, be used to tune the total x^ 3 ) response of 
the system. 

In our analysis, we have to consider boundary conditions which lead to uncountably many modes. In addition, in 
PCs the resonant set of arguments cannot be unambiguously separated from the nonresonant ones since the nontrivial 
resonances have to be determined numerically and the influence of the refractive index contrast on their positions 
cannot be quantified (see Sec. IIII[) . Similar to the cases of finite degrees of freedom with purely imaginary spectra of 
incommensurable eigenvalues (depending on whether or not un-removable resonances remain, the so-called Birkhoff 
of Birkhoff-Gustavson normal forms) we utilize - for our case of uncountably many degrees of freedom - normal forms 
as a formal tool to simplify the Hamiltonian and, hence, the equations of motion, but leave aside certain questions 
regarding the convergence of the transforms. 

A. Normal form transformation - Lie series 

In view of the discussion pertaining the efficiency of wave interaction processes in conservative systems (see Sec. 
IIIBI) . we assume in the following that the initial conditions, i.e. the initial distribution of modes dfc(O), a£(Q) in the PC 
with quadratic nonlinearity, are such that the resonance condition of lowest order that can be satisfied is a four-wave 
mixing process of the form 

L0 kl + u k2 - L0 k3 - u ki = 0. (45) 

All other resonance conditions that may hold are assumed to be of order 5 or higher. This four-wave resonance 
condition, Eq. (|45|) . is trivial in the sense that it can be fulfilled for any kind of dispersion relation and all possible 
initial conditions, even if the spectrum does not include the frequency w = 0. Consequently, the system Hamiltonian 
H = H2 + i?3 + H4 + ■ ■ ■ comprises contributions in H3 and H4 that lead to oscillatory source terms in the equation 
of motions for the time-dependent amplitudes a k {t), a\{f), as long as these modes are bound away from the resonant 
subsets. To simplify the structure of these equations, we search for a coordinate transformation that eliminates all 
terms in the Hamiltonian that correspond to resonance conditions which cannot be fulfilled. In constructing such a 
transformation, we follow the ideas outlined in [lH [T^. [l7j. 

In order to split the Hamiltonian into resonant and nonresonant parts, it is useful to introduce the linear mapping 

ad F (G) = {F,G}, (46) 

where {F, G} is the Poisson bracket of two functions F, G in phase space. In quantum optics, the operator adp acting 
on any function G is known as the Liouville operator. As an illustration, let us consider the Hamiltonian presented in 
Sec. HIl It can be easily verified that for ' the following relation holds: u kl + uj k2 +uj k3 = V(fci, k^, ks) € R 3 =£■ 

adH 2 {H^) — 0. More generally, adu 2 splits the Hamiltonian into resonant H™ s and nonresonant terms of order 
m, where H™ s e ker(ad H2 ) = {F m : ad H2 (F m ) = 0}. 

Now, our goal is to eliminate all nonresonant terms from the Hamiltonian and to replace the initial system containing 
resonant and nonresonant terms with a simpler system. In particular, we want to find a canonical (symplectic) 
transformation d k — F(b k ), such that in the new variables b k of the transformed Hamiltonian H'(b) takes on the form 

H'(b) = H 2 (b) + H'r s (b) + R'i(b), (47) 

where the residual R'^(b) contains only processes of fifth or higher order in the amplitudes b k . Such reductions can 
be performed by means of the Lie series transformations, which in the nonlinear optics community received attention 
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through the concept of the guiding center soliton in amplified fiber lines [10]. We define the Lie series exp(adx) 
induced by the generating functional K as 

exp(ad K )(H) = V —ad™{H). (48) 
^— ' to! 

m=0 

Here, X is a minimally third-order functional in the dynamic variables. For the purposes of the present work, only a 
few of the remarkable properties of these series will be used. Among them are [1a [191]: 

exp(ad K ){F,G} = {exp(ad K )F,exp(ad K )G} (49) 
(exp(adx)) 1 = exp(ad^x), (50) 
exp(adK)F(u) = F(exp(ad K )u), (51) 

where F, G are functions in phase space and u are the phase space variables. The first equality establishes the property 
of being a symplectic transformation for arbitrary generating functional K and can be verified by using the Jacobi 
identity for the Poisson bracket. The second equality makes a connection to the inverse Lie series transformation. 
Finally, the last equation allows us to interchange the order of mapping. 

The explicit form of the polynomial K can be determined by additional requirements on the normalization of the 
Hamiltonian. Writing K — K 3 + K4, where K 3 ,K±, respectively, are - yet to be determined - functionals of order 3 
and 4 in the dynamical variables, and expanding the exponential yields 

exp{ad K )H = exp(ad K ) [H 2 + H 3 + H{ cs + H% T + R 5 ] 
= H 2 + H 3 + H T 4 es + Hl T + R 5 + {K, H 2 } 

+{K,H 3 } + ±{K,{K,H 2 }} + R' 5 . (52) 

This suggests that the choice 



K 4 = ad H \ ( [H 4 + lad H3 (K 3 )] m ) (54) 



K 3 - ad H \{H 3 ), (53) 
1 

2' 

realizes the desired simplification, i.e., we obtain 

exp(ad K )H = H 2 + H{ es + B%. (55) 

In Eq. (|54p . ]-] nr indicates that we only consider the nonresonant terms contained within the bracketed expression. 
Alternatively, we can apply a transformation from the original set of coordinates a^, a* k to a new set bk, b* k according 

to 

a k = exp(ad K{b ))b k , (56) 

where K is the functional defined in Eq. (|54[) and we have - and will continue to do so in the remainder of this paper 
- dropped the bars over the coordinates. More explicitly, we choose 

K 3 = [ dk 123 A k % 2k3 b kl b k2 b k3 S(G -ki-k 2 - ha) (57) 

•* G 

+ [ dk 123 A k 2) ik2k3 b* ki b k2 b k3 S(G + h - k 2 - k 3 ) + c.c, (58) 



G 



as well as 



Ki = dk 123A B { k ^ k2kaki b kl b k2 b k3 b ki Y S ( G -h-k 2 -k 3 - ki) (59) 

J G 

+ / dk 123A B { k % 2kaki b* ki b k2 b k3 b ki ^5{G + k 1 -k 2 -k 3 - fc 4 ) + c.c. 
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One can note that since ({F, G} a , a *)* — {F* ,G*} a , a * , the requirement K = K* follows immediately from the 
definition of the coordinate transformation. The unknown coefficients A^' in K% are determined from the nonresonance 
condition, Eq. ([54]). and are readily found to be 

A ^k s = i + T k3 +UJ > (6°) 
Wfcj + Wfe 2 + ^k 3 

a( 2 ) - A k t k 2 k 3 

Equation (|54[) also fixes the coefficients B^ % \ which can be expressed in terms of and AW, i — 1,2. However, 
their explicit form is not required if we are only interested in terms at fourth order to which K4 does not contribute. 
Finally, we note that the expressions for become singular if the resonance condition is fulfilled. 
In these new coordinates, the transformed Hamiltonian reads 

H'(b) = H 2 {b)+H'r s {b)+R'l{b) 
= J dku k b k b* k 

+ / dk 12 3iS klk2k;ski b kl b k2 b k H b k ki ^25{G + fei + k 2 - k 3 - k 4 ) 

•* G 

+B!£(b). (62) 

With the help of the corresponding Poisson brackets 

{b k ,b* k ,} = iS(k-k'), (63) 
{b k ,b k ,} = 0, (64) 

we finally obtain the transformed equations of motion 

—b k = tuj k b k 
dt 

+i / dk 123 (S klk2k3k + S klk2kk3 )b kl b k2 b* k;i y^ y S(G + kg + fc 2 - fc 3 - k) 

•* G 

+R±(b), (65) 

where R^b) is a residual containing higher-order terms. 

Equation (|65[) represents the generalization of the Zakharov equation to the case of periodic media. This equation, 
which Zakharov reduced in his original article [f| to the Nonlinear Schrodinger equation (NLSE), has been the subject 
of extensive studies for waves in homogeneous systems and has been tested experimentally for the case of surface 
water waves in channels 



B. Effective fourth-order coefficient 

Using the coefficients A^ l \ we now proceed to calculate {K(b), H 3 (b)} and its contribution to the transformed 
Hamiltonian, which contains the renormalized fourth-order matrix element. Indeed, the effective fourth-order element 
S kl k 2 k 3 k 4 consists of the term originating from the third-order nonlinear susceptibility of the system and three 
additional terms i = 1,2,3 that arise from nonresonant \^ processes 

S kl k 2 k3ki = ^kik 2 k 3 k 4 + 2 {^kX 2 k 3 k 4 + ^kik 2 k 3 ki + ^kik 2 k 3 kij (^6) 
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where the explicit expressions for the are 

(v {1) 

iy(l) _ | {-G-k 1 -k 2 )k 1 k 2 (-G-k 1 -k 2 )k 3 k 4 f (-G-fc 3 -fc 4 )fc 3 fc 4 (-G-fc 3 -fc 4 )fcifc 2 

iV fcifc 2 fc 3 fc 4 ~ ~ J 2 / 



y-(l) -5/(1)* -j/fl)* 



•'(-G-fci-fca) + + W fc2 W(_ G _ fca _ fc4 ) + W fc3 + CJ fe4 

t/(2) t/(2)* T/ (2)* T/ (2) 

jy(2) _ \ "* ( r (G+fci+fc 2 )fcifc 2 V (G+fcl+fc 2 )fc 3 fc 4 _,_ ^(G+fc 3 +fc 4 )fc 3 fc 4 ^(G+fc 3 + fc 4 )fcifc 2 



k t k 2 k 3 k 4 



J k\k 2 k 3 k 




->(G+fci+fc 2 ) + ^fei + w fc 2 -^(G+fe 3 +fc 4 ) + Wfc 3 + U7fc 4 



T/ (2) T/ (2)* T/ (2)* „(2) 

(3) . I "fc 3 (-G+fc 3 -fei)fci fc 2 (-G+fc 3 -fci)fc 4 . %(-G+fci-fc 3 )fe 3 r fc 4 (-G+fci-fe 3 )fe 2 



-UJ k3 + ^(-G+fc 3 -fci) + w fci —&ki + w (-G+fei-fe 3 ) + <^fc 3 



In the above expressions, the denominators of N kik2k ki contain the resonance conditions. In particular, frequencies 
which depend on the difference between two wave vectors appear. Note that Sk 1 k 2 k 3 k i describes coupling of the four 
waves bki , and therefore, trivial resonances appear in the denominators for any possible distribution of wave vectors 
in the wave packets b ki . As a consequence, Sk lk2 k 3ki is ill defined, unless the singularity is removable. To analyze this 
issue in the next section, we further simplify the matrix elements by considering the simplest possible situation, i.e., 
the propagation of a monochromatic wave. 



C. Simple case - monochromatic wave 



In order to obtain further insight into the physics behind the coefficient, and to discuss the technical issue related 
to the appearance of resonances, we consider the simple situation when initially only a single monochromatic wave is 
excited in the transformed space 



b k = BS(k-k ), 

where B — const. £ C and fco ^ 0. In this case, the expression for the effective fourth-order element reduces to 



(67) 



Sk k k fco 



(3) 



fco fco fco fco 



V, 



(1) 



(— G— 2ko)koko 



E 

G 



(2) 



(G+2fc )fc fco 



V, 



-^G+2fc + 2Wfe 



-4 



(2) 



fco(-G)fc 



G 



(68) 



At this point, the following comments are in order: 

(i) The terms arising from the nonrcsonant interaction processes include a summation over all reciprocal lattice 



vectors G. Since the coefficients V, 



(i) 



fci fc 2 fc 3 



and frequencies 10 ki are evaluated in the extended zone scheme, a 
summation over reciprocal lattice vectors G is equivalent to a summation over all band indices n. 

(ii) For G = 0, the last term is potentially undefined since luq = 0. Fortunately, this singularity is removable. 
Indeed, as discussed in Sec. [Til the Bloch functions are normalized as D k cx y/Uk. Therefore, we obtain 



lim 



(2) 



fcoM(fc'o+p) 



lim 



<C(fc o ,0), 



(69) 



with some constant C(fc,0) for which < C(k,0) < 00 V(fc) € 



(iii) The denominator of the terms in the second line of Eq. (|68l) include the squared modulus of the initial three- 
wave matrix element. From Eqs. (|31[) and (|32[) it is clear that the nonlinear susceptibility enters with a power 
of 2 into the numerators. We would like to note that the resonance conditions appearing in the denominators of 
Ni correspond precisely to the three-wave interactions that have been eliminated from the Hamiltonian. As a 
result of this elimination, these processes give effective contributions to the interactions arising from nonlinear 



coupling via x^- Considering the contribution of N^* J k k k , the phase mismatch of second harmonic generation 



r(2) 



appears in the denominator, whereas for N-. 



(3) 

fc'o fco fco fco 



, the sum of frequencies accounts for the interaction of two 
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waves with the same frequency via the static field luq or via waves at symmetry points log- In the field of 
nonlinear optics of dispersive homogeneous media, effective higher-order nonlinear contributions due to phase 
mismatch are generally referred to as cascaded processes. These effects are known to effectively modify the 
nonlinear third-order response of homogeneous materials with noncentrosymmetric nonlinear susceptibilities 
[20I ] . The phase mismatch in these materials results from intrinsic material related dispersion, which - due 
to Kramers-Kronig relations - always appears together with absorption. Although some progress has recently 
been made on the Hamiltonian description of the optical properties of nonlinear dispersive materials j2l| , the 
Lie series transformation does not directly apply to this class, because the detailed knowledge of interaction 
coefficient relies on a modal description. The key aspect about PCs which enable the formal construction of 
a normal form is that nontrivial dispersion relations result from the formation of a photonic band structure. 
Consequently, if the validity of the equation could be justified for finite times with an appropriate error bound, 
the result suggest that it is possible to achieve cascading without using strongly dispersive and, therefore, lossy 
absorbing materials. In addition, the judicious engineering of PCs would allow to tailor the effective third-order 

nonlinear response of these systems through both the dispersion relation u>k and the matrix elements VPl k 

(2) 

and Vfc h k . Taken together, these two aspects enlarge the possibilities for the realization of optical materials 
with tailored nonlinear properties. 

Although some results on the validity for finite times of envelope equations that have been derived under nonresonance 
assumptions in the presence of resonances are known [22[, these turn out to be inapplicable to the present case, as 
outlined in the following section. 

D. Influence of resonances 

As a consequence of the strong asymptotic resonance 

r(k, to) = u>(k) — uj[k — to) — oj(to) — > (70) 

for fc, to 0, the system truncated after the third order terms cannot be expected to be well posed. In detail, we have 
ujk ~ k plus some correction terms coming from the periodic coefficients. It is well known [23[ that the corrections 
decay asymptotically with the regularity of the associated coefficients. If the coefficients are n-times differentiable, 
then for suitably chosen k and to, we obtain 

r(k,m) < C (\k\~ n + \k - m\-' n + \m\- n ) . (71) 

Since r(k, to) appears in the denominator in the normal form transformation, there is a substantial loss of regularity 
associated with the normal form transformation. Thus, the less regular the coefficients, the better are the conditions 
for well posed truncated third-order systems. 

The same problems occur if the truncated third-order system is justified with some approximation result similar to 
that of Ref. [2ij . In this work, the NLSE has been justified in the above sense for semi linear wave equations with 
periodic coefficients. An additional difficulty occurs in our situation due to the quasilinearity of the problem that 
always leads to some loss of regularity in normal form transformations. 

However, on physical grounds (see Sec. IIIip we expect that for suitably chosen initial conditions, the dynamics of 
the system is dominated by nonresonant processes for not too long times. Quantifying this conjecture by methods of 
asymptotic analysis is a challenging mathematical problem that has not been addressed yet. 

V. CONCLUSIONS 

In this article, we have extended a recently developed Hamiltonian formalism to the analysis of nonresonant wave 
interaction in one-dimensional PCs with quadratic nonlinearity. In particular, we have shown that an appropriate 
normal form transformation leads to a generalization of Zakharov's equation for periodic media. For the simplest case 
of monochromatic wave, the resulting effective coupling coefficient in the transformed equations of motion is the result 
of cascading second-order nonlinear processes. Moreover, we derived explicit expressions for the emerging effective 
nonlinear coupling terms. This may allow one to tailor the effective nonlinear response of the composite materials 
through judiciously engineering of its linear properties and could, in turn, suggest possibilities for the realization 
of optical materials with customized nonlinear properties. The approach presented here can also be applied to 
nonresonant higher-order interaction processes. Owing to the linear dispersion relation in the long-wavelength limit 
certain mathematical difficulties arise that limit the applicability of the method to finite times. However, we would 
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like to emphasize that in experiments, various loss mechanisms or the finite sample size set maximal time scales so that 
this limitation may be less severe. In addition, full numerical simulations can provide further insight into the range of 
validity of our theory. A systematic study of the Zakharov equation's validity represent a challenging mathematical 
task and further research is indicated. 
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